Van der Waals CoefRcients of Atoms and Molecules from a Simple Approximation for 

the Polarizability 
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A simple and computationally efficient scheme to calculate approximate imaginary-frequency de- 
pendent polarizability, hence asymptotic van der Waals coefficient, within density functional theory 
is proposed. The dynamical dipolar polarizabilities of atoms and molecules are calculated starting 
from the Thomas-Fermi-von Weizsacker (TFvW) approximation for the independent-electron ki- 
netic energy functional. The van der Waals coefficients for a number of closed-shell ions and a few 
molecules are hence calculated and compared with available values obtained by fully first-principles 
calculations. The success in these test cases shows the potential of the proposed TFvW approximate 
response function in capturing the essence of long range correlations and may give useful information 
for constructing a functional which naturaUy includes van der Waals interactions. 

PACS numbers: 



I. INTRODUCTION 

Density functional theory (DFT)i has become a stan- 
dard tool for electronic structure calculations of atoms, 
molecules and materials thanks to the widely used lo- 
cal density approximation (LDA)^ and generalized gra- 
dient ones (GGA)'^ to the exchange-correlation function- 
al in the Kohn-Sham (KS) formulation. These approxi- 
mations describe well many properties-such as cohesion, 
bonds, structures, vibrations, etc.- of densely packed 
molecules and materials when covalent, metallic or hy- 
drogen bonds are involved. However, they fail when 
applied to systems having regions of small overlapping 
density where long range correlation effects, which are 
not treated correctly in LDA and GGA, are important. 
These kinds of systems are frequently met in nature, 
e.g. in biomolecules, as well as in physical and chem- 
ical processes, such as, for instance, in molecular ad- 
sorption on surfaces, chemical reactions, etc. The sim- 
plest paradigmatic examples of these weakly bonded sys- 
tems are dimers of noble gas atoms and/or closed-shell 
molecules and several works in the literature exists where 
the van der Waals (vdW) Cq coefficients are calculated 
within DFT at different levels of sophistication, ranging 
from a crude treatment of response functions in term of 
electronic densities^ to full calculations within the time- 
dependent DFT framework.^ 

During the past decade, many attempts have been 
done to improve performances of DFT calculations in 
these systems j^iiiSii^ by including explicitly in the ex- 
change and correlation functional non-local correlations, 
that are missing both in LDA and in GGA. The general 
staring point in all these approaches is the so-called Adia- 
batic Connection Fluctuation-Dissipation (ACFD}i2, ex- 
pression for the exchange-correlation (xc) energy, where 
an exact formula is used, giving this quantity through dy- 
namical response functions of all fictitious systems which 
connect the non-interacting KS system with the real 



many-body interacting one. These calculations are how- 
ever extremely demanding for practical applications and 
have been performed only for a limited number of sys- 
tems so far—. Recently, an approximate scheme has been 
proposed where the ACFD was used as starting point for 
further simplifications and which has demonstrated to be 
able to account for van der Waals interaction with some 
success for a number of cases.— 

In our opinion, more calculations and further devel- 
opments need to be pursued in order to explore the po- 
tential offered by an approximate treatment of ACFD for 
the calculation of accurate exchange-correlation energies. 
Our aim in this work is to assess the possibility of ac- 
curately approximate the frequency dependent response 
function of a system starting from the Thomas-Fermi 
and von Weizsacker (TFvW) approximation for the non- 
interacting kinetic energy functional. Only the dipolar 
response function responsible for the asymptotic disper- 
sion interaction between non overlapping fragments will 
be considered here. The encouraging results obtained for 
the Cg coefficients will give support to the possibility of 
using similar procedures in the extremely demanding cal- 
culation of accurate correlation energy based on ACFD. 

The use of approximate kinetic energy functionals to 
evaluate dynamical response properties of materials is 
not a new ideaiS but is constantly appealing one and, for 
instance, Barejee and Harbolai^ have recently calculated 
approximate vdW coefficients for large alkali metal clus- 
ter using an TFvW approximation in conjunction with 
an hydrodynamical approach to the collective excitation 
of the system. 

In this study we follow a different approach to the 
same problem, where, starting from Thomas-Fermi and 
von Weizsacker approximation for kinetic energy func- 
tional, we develop a simple and computationally fast 
procedure to calculate imaginary-frequency-dependent 
polarizability-hence the van der Waals interaction in the 
asymptotic region-by Density Functional Perturbation 
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Theory (DFPT)^ii^ Our method contains the follow- 
ing ingredients: (i) the ground-state (GS) electronic den- 
sity, n(r), which is accurately computed from the stan- 
dard KS procedure within LDA or GGA; (ii) a single aux- 
iliary wave function, ip{r), corresponding to the GS den- 
sity, and the auxiliary Hamiltonian admitting 'p{r) as its 
GS eigenfunction in the TFvW approximation; (iii) the 
so-called modified Sternheimer equation which is used 
to calculate the imaginary-frequency polarizability of the 
system. This procedure can be considered as a simplifica- 
tion of the one proposed by Mahan^^ and others^^'^^ more 
than twenty years ago for the exact- within LDA/GGA 
DFT-calculation of atomic and molecular polarizabili- 
ties. The original scheme becomes very demanding for 
systems having a large number of electrons while our ap- 
proximate treatment, dealing only with one single aux- 
iliary wavefunction, can deal efficiently with systems of 
any number of electrons. 

This scheme could be a good starting point, we be- 
lieve, for developing an approximation of the xc-energy 
functional which is simple enough to allow broad appli- 
cations but still capture the essence of van der Waals 
energies. Moreover, it could also be useful for calculat- 
ing vdW coefHcients of fragments which are needed in 
some semi-empirical treatments of vdW energy in DFT 
such as the damped-dispersion force method . ^^'^^ 



II. THEORY 

A. Calculation of van der Waals coefficient within 
DFPT 

From the theoretical point of view, van der Waals co- 
efficients are most conveniently determined through the 
frequency-dependent polarizability, a{iu), of the system 
thanks to the relation 



du ai{iu)a2(iu). 



(1) 



In the framework of density functional theory, the static 
polarizability, a{0) can be obtained from the density re- 
sponse of the system under a uniform external electric 
field (along the z-axis) via 



a(0) ^ 2 



zAn(r) 
E 



dr. 



(2) 



An efficient computational technique for calculation of 
density response that avoids the computation of unoccu- 
pied states of the independent-particle Hamiltonian has 
been proposed by Sternheimer^^ more than fifty years ago 
for atomic polarizability calculations. The method was 
later modified by MahaniS, to include a self-consistent 
treatment of the electrons for the calculation of atomic 
polarizabilities within density-functional theory in LDA. 
The linear density-response. An, to an external perturba- 
tion, AVext, is determined by the following self-consistent 



set of linear equations. 

An = 2^7^e[V'*AV'i] 



(3) 

AVks = AV^^t + AVh + Av^c, (4) 
Ai,,^-\AVKS- Ae,\'4,,, (5) 



V2 

— ^ + Vks - e^ 



where the sums run on the set of occupied orbitals only. 

A generalization of the method for calculation of po- 
larizability at a finite imaginary frequency (see below) 
was also made by Mahan and was used to calculate 
the van der Waals coefficients for a number of atomic 
systems.—' The method has later been adapted and suc- 
cessfully applied to the case of extended periodic systems 
where it is better known as Density Functional Perturba- 
tion Theory. This approach is exact within LDA DFT, 
however, computationally demanding since the computa- 
tional cost grows as the third power of the system size. 



B. Calculation of van der Waals coefficient using 
TFvW functional 



One possible way to reduce the computational cost in- 
volved in the DFPT calculation of static and dynamical 
response functions is to approximate the non interacting 
kinetic energy with some orbital free functional, as in the 
TFvW scheme. 

Let us assume that the GS density distribution, n(r), 
of an atomic or molecular system has been computed ac- 
curately within the KS scheme, employing, for instance, 
LDA or GGA as xc-functional, and let us introduce the 
auxiliary wave function (/^(r), which is normalized to unit 
and simply related to the electron density as 



n(r) = iV|(p(r)|^ 



(6) 



where iV is the number of electrons in the system. The 
crucial point in our scheme is the assumption that the 
response of the system around its GS density can be ap- 
proximated by using the TFvW functional for the kinetic 
energy. In this approximation the total energy functional 
in term of ^piy) reads 



Exciv] + Eext[^] - N^i 



\ip\^dv-l 



(7) 



where a = ^(Stt^)^/'^ and /i is the Lagrange multiplier 
used to enforce normalization (Rydberg atomic units are 
used throughout the text). The choice of this approx- 
imation is inspired by the fact that the vW correction 
term gives the exact kinetic energy in regions where only 
one wave function is relevant, typically the asymptotic 
region of atoms or molecules. Moreover, the dominant 
contributions to the polarizability come from the loosely 
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bound electrons in this region, that are thus expected to 
be captured in this approximation. 

The corresponding Euler equation of Eq. ([7]) deter- 
mines fir) 



2 



ip{r) = 0, (8) 



where kp{r) — (37r^n(r))3 is the local Fermi wave- vector. 
If Vexti^) is given, an approximate GS density in TFvW 
approximation can be obtained by solving this equation. 
Here we invert the reasoning and the density, hence ip(j), 
is assumed to be given and an auxiliary effective poten- 
tial, denoted by 14// (r), is constructed such that the 
corresponding Hamiltonian admits (p{r) as its GS eigen- 
function. It can be formally written as 



T/ ( A 1 VV(r) 



(9) 



and Ve//(r) can be found from this equation once (f{r) 
is known, even if some care must be taken in the asymp- 
totic region where the density, and hence ip(j), vanishes 
exponentially. 

The linear density-response to an external perturba- 
tion is determined by the following self-consistent set of 
equations. 



An = 2NTZe[ip* A(p], 



AVeff = AVext + AVh + Av^c + '-T^An 



V2 

— ^ + K// - M 



Alp ■ 



(10) 

(11) 

[AV^ff - A^]^, (12) 



/r2 

in 



where A/x = {Lp\AVef f\f) , and Aif, which satisfies the or- 
thogonal condition {lp\Alp) = 0, is the first order change 
of the wave function brought about by the external field. 
Once this set of equations is solved, static polarizability 
is calculated from Eq. In order to find the polar- 
izability at finite imaginary frequency, a{iu), the above 
procedure is slightly modified by adding a frequency term 
iu to /Lt, making it become a complex quantity /i -|- iu}'' 
From the imaginary- frequency polarizability, a{iu)^ the 
vdW coefficient Cg can be immediately calculated from 
Eq. Q. 
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FIG. 1: Dynamic polarizabilities of several arbitrarily chosen 
ions calculated by our method compared to corresponding 
quantities in Ref. [l7|. 



be obtained directly from Eq. For the latter, the situ- 
ation is more complicated because all quantities are cal- 
culated within the plane-wave pseudopotential method 
and are expanded in Fourier components up to a given 
kinetic-energy cut-off which makes the inversion needed 
in Eq. [H] numerically difficult. We have overcome this 
difficulty by an iterative optimization process inspired by 
the method described in Ref. for the determination 
of the optimized effective potential starting from an ac- 
curate density: assuming that at a given iteration, i*'*, 
the approximate effective potential is V^jf{v) a residual 
quantity ^^(r) is defined by 



2 



+ ^e//W-A* 



^(r). (13) 



This quantity vanishes everywhere only if the Hamilto- 
nian corresponding to this potential admits <p(r) as its 
GS eigenfunction. As long as this is not the case the 
potential is updated as 



K7/(r) = K//(r)+a^^(r)+/3, 



(14) 



where a and [3 are chosen in such a way that the norm of 
the new residual ||5'*'^^|| = J dr [S''+^(r)]^ is minimized. 
The process is terminated when the integrated charge 
difference Sn'' — ^ f dr |n*(r) — n(r)| is less than a given 
threshold, practically chosen to be of the order of lO^'^ . 



C. Construction of the effective potential 

We have chosen to apply our method to calculate dy- 
namic polarizabilities and vdW coefficients for two kinds 
of systems, namely spherically symmetric ions and some 
simple molecules to demonstrate the efficiency of the 
method. By exploiting the symmetry properties of the 
former systems, the GS density n(r) can be calculated 
very accurately by integrating the radial KS equations 
on a logarithmic grid with the highly accurate Numerov's 
algorithm. In this case the effective potential T4//(r) can 



III. RESULTS 
A. Atomic systems 

Numerical results for closed-shell ions show a good 
agreement with those obtained by Mahan with an ac- 
curate but more expensive calculation for a wide range 
of atomic number and frequency. This can be seen 
in Fig. 1 where frequency-dependent polarizabilities of 
some closed-shell ions obtained by the two methods are 
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TABLE L Ce values for dimers (Rydberg atomic units). 
Present: our results. Mahan: calculations from a very sim- 
ilar scheme in Ref. il7i . Unified: calculations from Ref. [271 
using self-consistent electrodynamics. Reference: values cited 
as reference in Ref. [27l . 



Dimer 


Present 


Mahan 


Unified 


Reference 


He 


2.64 


3.64 


2.58 


2.92 


Ne 


15.44 


13.96 


15.0 


13.8 


Ar 


133 


132.2 


143 


134 


Kr 


266 


261.4 


291 


266 


Xe 


600 




663 


597 



compared. In Table |T] calculated values of vdW coeffi- 
cients for a number of pairs of rare gases are given. Our 
results are all in the range and at least as accurate as 
those reported in Ref. i27, when compared to the refer- 
ence ones. 

In Fig. 2 our calculated Cg values both for homonu- 
clear and mixed pairs of 14 ions are plotted against those 
reported in Ref. The good agreement between the 
results obtained by the two methods are indicated by the 
narrow spread of the points around the diagonal. Quan- 
titatively, the difference never exceeds 25% and in most 
cases is less than 10%. 




Our calculated C, coefficients 

o 

FIG. 2: Ce values of all possible pairs of 14 ions calculated 
by our method plotted against corresponding values shown in 
Ref. [13 



B. Importance of using a good charge density 

The calculation of van der Waals coefficients presented 
above is from response functions in TFvW approxima- 



tion calculated around an accurate charge density. For 
this purpose, we need to introduce an effective poten- 
tial admitting the square root of the charge density as 
its ground state wave function as already discussed in 
Sec. II B. One may wonder if this is indeed necessary 
since it seems natural to calculate the charge density 
also in TFvW approximation and then use it as input 
for the calculation of Cg coefficient. Moreover, doing 
calculation in this way makes the construction of the ef- 
fective potential needed in our calculation unnecessary 
because it is determined in the self-consistent solution of 
Eq. ([5]). We have tried this option for the case of no- 
ble gas atoms and the results is disastrous. For example, 
Cg coefficient of He changes from 2.1 a.u. when com- 
puted with accurate LDA-DFT charge density to 15.1, 
and 227.2 (a.u.) when calculated from the charge densi- 
ties obtained from the solution of the Hartree equation 
or from solving self-consistcntly the TFvW approxima- 
tion. Reducing the weight of the gradient corrected term 
in the functional to 1/5, an empirical value often used 
in the literatur e — t^-, still gives a very poor result (36.4 
a.u.). This behaviour is not totally unexpected since it 
is well-known^^ that also the TFvW kinetic energy itself 
- the quantity on which the approximate response func- 
tion is based - while giving accurate estimates when ap- 
plied to accurate charge density behaves poorly if treated 
self-consistently. This result indicates the importance of 
calculating the response functions with accurate charge 
densities and shows that our approach, though not being 
a self-consistent procedure, is the correct way to calculate 
vdW coefficient using TFvW approximation. 



C. Role of core electrons 

For systems without spherical symmetry, the KS equa- 
tions as well as the modified Sternheimer ones are no 
longer radial-angular separable, and a general method 
to calculate GS electronic structures must be used to 
solve them. The scheme described above have been im- 
plemented in the PWscf plane-wave pseudopotential code 
which is part of the Quantum ESPRESSO distribution 
Although in the pseudopotential approach the core- 
electron density could be easily included in the defini- 
tion of the auxiliary function, i^(r), it is not convenient 
because of the higher computational cost due to larger 
kinetic-energy cutoff that would be required in the cal- 
culation of the GS density n{r). It is therefore worthwhile 
to estimate the contribution, expected to be small, of the 
core electrons to the polarizability in this scheme. To this 
end, we compared the results obtained in atomic calcu- 
lations where (p{r) was computed from the total charge 
density or from the valence-only charge density or, as it 
is done in the non- linear core correction (nice), from the 
valence charge density plus a smoothed core charge. Here 
we show the results for Beryllium and Argon atoms for 
which the effects of the core charge on the total polariza- 
tion are expected to be considerably different. It could 
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FIG. 3: Static density responses of Be (top) and Ar (bottom) 
calculated with diflterent treatment of the core charge (see the 
insets): totally included (ae), completely neglected (ve), and 
partially included (nice). 



be not negligible for the former due to non-tightly bound 
core electrons to a small charge nucleus, while it should 
be very small for the later whose core electrons are more 
tightly bound. It turns out that in both cases the contri- 
bution of the core charges to Cg values is very small. In 
Fig. [3] static response densities calculated by using differ- 
ent densities in the core region (see the insets) are plot- 
ted. This figure shows two opposite roles of core electrons 
in determining the polarizability. On the one hand, they 
contribute to the density response in the core region, thus 
making the total polarizability larger. They prevent, on 
the other hand, penetration of valence electrons into the 
core region when an electric field is applied, therefore re- 
ducing the total polarizability. Although these contribu- 
tions do not cancel each other exactly, they still make the 
effects of the core charge not very important as indicated 
by the values of static polarizabilities and Cg coefficients 
reported in Table [Til where the difference between totally 
included and completely neglected core charge is just a 
few percent. Therefore, in the calculations for molecules 
presented below, a little accuracy has been sacrificed by 
using only the density of valence electrons in order to 
reduce the number of plane-wave needed to describe the 



TABLE IL Static polarizabilities a(0) and vdW coefficients 
Ce (Ry atomic units) of Be and Ar calculated with differ- 
ent densities: all-electron (ae), valence with some core-charge 
(nice), and only valence (ve) density. 







a(0) 






Co 






ae 


nice 


ve 


ae 


nice 


ve 


Be 


33.07 


33.05 


33.54 


194.10 


194.14 


190.58 


Ar 


10.93 


10.92 


10.96 


66.80 


66.85 


63.80 












-- Full 






TF+vW 















Frequency (a.u.) 

FIG. 4: Imaginary-frequency dependent polarizabilities of 
methane (the inset) and benzene molecules calculated by 
TFvW method (solid red curves) compared to results of full 
calculation (dashed black curves). 



auxiliary wavefunction of the system. 



D. Molecular systems 

To exemplify the general scheme, we have calculated 
dynamic polarizabilities and vdW coefficients of methane 
and benzene, two molecules with different nature of 
chemical bonds and geometric structures. The KS equa- 
tions for each isolated molecule were solved using peri- 
odic boundary conditions in a simple cubic simulation cell 
with side length of 12 and 10 A and kinetic-energy cutoffs 
of 80 and 60 Ry, respectively. Simple LDA exchange- 
correlation functional with norm-conserving pseudopo- 
tentials were used to obtain the GS charge density of the 
isolated molecules. 

Fig. [3]shows the imaginary- frequency dynamic polariz- 
abilities of methane and benzene molecules calculated in 
our scheme, compared with the result of the full calcula- 
tion which has also been implemented in the PWscf plane- 
wave pseudopotential code. For methane molecule, the 
result of our simplified calculation compares excellently 
with the one of the more accurate method. Although 
this is not the case for benzene molecule, nevertheless 
the difference between the two calculations is still rather 
moderate as expected. 

From the dynamical polarizability the Cg coefficients 
can be obtained and the approximate calculation agrees 
very well with the full calculation for methane (264 a.u. 
to be compared with 271 a.u. obtained in the full calcu- 
lation) while for benzene TFvW approximation overes- 
timate by about 40% the full calculation (4.9 xlO'^ a.u. 
w.r.t. 3.6 xlO^ a.u.). 
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IV. DISCUSSION AND CONCLUSION 

In this study we have shown that a simple approxi- 
mation for the (dynamical) polarizability of atoms and 
molecules can be obtained starting from the TFvW ap- 
proximation to the independent electron kinetic energy 
functional. The success of the present approximate 
method in the calculation of frequency dependent po- 
larizability and vdW coefficients of atomic and molecular 
systems, together with its computational efficiency, espe- 
cially for large systems, suggests that it may be a useful 
tool to explore the behavior of systems where non-local 
long-range correlations are important. 

An algorithmic development of time-dependent den- 
sity functional theory in real time has recently been pro- 
posed by Marques and co-workers2^ that allows for an 
efficient evaluation of van der Waals coefficients from 
the full Kohn-Sham response function. The computa- 
tional cost for the calculation of Cg coefficients in this 
approach grows quadratically with the system size. In 



our approximate scheme instead the computational cost 
grows only linearly with the size of the system. This is 
because only one auxiliary wave function is needed no 
matter how many electrons are present in the system. 
Even for the small atomic and molecular systems con- 
sidered in this work the computational time required by 
the simplified calculation is at least one order of mag- 
nitude lower than that of the full calculation and this 
is expected to become increasingly more convenient for 
larger systems. Only a systematic comparison of the re- 
sults obtained for more systems within our approximate 
scheme as well as within other approaches will allow to 
better assess the computational and physical merits and 
limitations of the various approaches. 

One of the authors (HVN) would like to thank the 
Abdus Salam ICTP, where the initial part of this work 
was done, for financial support in the framework of 
ICTP/SISSA Joint Master's Degree Programme (2003- 
2005) and the hospitality at the centre. 



^ P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). 
2 W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965). 
^ J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett. 
77, 3865 (1996). 

Y. Andersson, D. C. Langreth, and B. I. Lundqvist, Phys. 

Rev. Lett. 76, 102 (1996). 
^ S.J. A. van Gisbergen, J.G. Snijders, E.J. Baerends, J. 

Chem. Phys. 103, 9347 (1995); V.P Osinga et ai, J. Chem. 

Phys. 106, 5091 (1996). 
® W. Kohn, Y. Meir, and D.E. Makarov, Phys. Rev. Lett. 

80, 4153 (1998). 

K. Rapcewicz, and N.W.Ashcroft, Phys. Rev. B 44, 4032 
(1991). 

^ J.F. Dobson, and J. Wang, Phys. Rev. B 62, 10038 (2000); 

J.F. Dobson, J. Wang, and T. Gould, Phys. Rev. B 66, 

R081108 (2002); 
^ H. Rydberg et al, Phys. Rev. Lett. 91, 126402 (2003); 

M. Dion et al, ibid. 92, 246401 (2004); ibid. 95, 109902 

(2005) ; 

^° D.C. Langreth, J. P. Perdew, Solid State Commun. 17, 

1425 (1975); D.C. Langreth, J.P Perdew, Phys. Rev. B 

15, 2884 (1977). 
" F. Furche, Phys. Rev. B 64, 195120 (2001); M. Fuchs and 

X. Gonze, Phys. Rev. B 65, 235109 (2002); A. Marini, P. 

Garcia-Gonzalez, A. Rubio, Phys. Rev. Lett. 96, 136404 

(2006) . 

^2 E. Zaremba, and H.C. Tso, Phys. Rev. B 49, 8147 (1994). 
" A. Banerjee, and M.K. Harbola, J. Chem. Phys. 113, 5614 

(2000); A. Banerjee, and M.K. Harbola, J. Chem. Phys. 

117, 7845 (2002); A. Banerjee, and M.K. Harbola, Pra- 



mana, J. Phys. 66, 423 (2006); A. Banerjee, and M.K. Har- 
bola, Phys. Lett. A 372, 2881 (2008). 

" R. M. Sternheimer, Phys. Rev. 96, 951 (1954). 

^5 G. D. Mahan, Phys. Rev. A 22, 1780 (1980). 

^® S. Baroni, P. Giannozzi, and A. Testa, Phys. Rev. Lett. 
58, 1861 (1987); S. Baroni, S. de Gironcoh, P. Giannozzi, 
and A. Dal Corso, Rev. Mod. Phys. 73, 515 (2001). 

" G. D. Mahan, J. Chem. Phys. 76, 493 (1982). 

M.J. Stott, and E. Zaremba, Phys. Rev. A 21, 12 (1980). 
A. Zangwill, and P. Soven, Phys. Rev. A 21, 1561 (1980). 
H.V. Nguyen, and S. de Gironcoh unpublished. 
X. Wu, M.C. Vargas, S. Nayak, V. Lotrich, and G. Scoles 
J. Chem. Phys. 115, 8748 (2001). 

G. Murdachaew, S. de Gironcoh, and G. Scoles, J. Phys. 
Chem. A 112, 9993 (2008). 
2^ S. Kiimmel and J.P. Perdew, Phys. Rev. Lett. 90, 043004 
(2003). 

^'^ E. Hult, H. Rydberg, B.I. Lundqvist, and D.C. Langreth, 
Phys. Rev. B 59, 4708 (1999). 

W. Stich, E.K.U. Gross, P. Malzacher, R.M. Dreizler, Z. 
Phys. A 309, 5 (1982). 

R.M. Dreizler and E.K.U. Gross, "Density Functional The- 
ory", Springer- Verlag, Berlin, 1990. 
^'^ E. Hult, H. Rydberg, B.I. Lundqvist, and D.C. Langreth, 
Phys. Rev. B 59, 4708 (1999). 

S. Baroni et ai, http://www.quantuin-espresso.org/. See 
also http://www.pwscf.org/. 

M.A.L. Marques, A. Castro, G. Malloci, G. Mulas, and 
S. Botti, J. Chem. Phys. 127, 014107 (2007). 



